Analyse Biome

DUBii 2021

Chargements des données et librairies

#installation de package nécéssaire en local
#install.packages("BiocManager")
#BiocManager::install("phyloseq")
#install.packages("phyloseq")
#install.packages("remotes")
#remotes::install_github("mahendra-mariadassou/phyloseq-extended", ref = "dev")
#install.packages("gridExtra")


#Installation des librairies
library(tidyverse)
library(phyloseq)
library(phyloseq.extended)
library(ggplot2)    ## for graphs  necessary ???
library(ape)        ## for tree manipulation
library(vegan)      ## for community ecology analyses

library(scales)

library(grid)
library(gridExtra)  ## for design graph
library(knitr)
library(kableExtra)

#chemin <- "~/save/Milint/Resultat/"

load("BIOM.Rdata")
#load("Data_Biom.Rdata")

Le BIOM a été créé dans le script01. Il est cependant “brut”.

Etape de filtre sur l’abondance par échantillon et OTU

Somme de tous les reads et on applique un filtre de 10-5 afin de s’assurer que ce ne sont pas des artefacts Dans notre cas permet de filtrer les taxa avec des valeurs supérieures à 180 reads par OTU

filtre <- sum(taxa_sums(Biom))*1e-5

filter_biom <- prune_taxa(taxa_sums(Biom)>filtre,Biom)

save(Biom, filter_biom, file = "Data_Biom.Rdata")

Description du BIOM

Niveau d’abondance par échantillon

# abondance sur données complètes
sample_sums(Biom) %>% range()

[1]   5042 240312

# abondance sur données filtrées
sample_sums(filter_biom) %>% range()

[1]   5003 239715

# moyenne 
moy <- sum(taxa_sums(Biom))/846
moyf <- sum(taxa_sums(filter_biom))/846
mi <- min(taxa_sums(Biom))
ma <- max(taxa_sums(Biom))
mif <- min(taxa_sums(filter_biom))
maf <- max(taxa_sums(filter_biom))
med <- median(taxa_sums(Biom))
medf <- median(taxa_sums(filter_biom))

df <- data.frame(
  échantillon = c(846,846),
  mean = c(moy,moyf),
  median = c(med,medf),
  min = c(mi,mif),
  max = c(ma, maf)
)
rownames(df) <- c("Biom","Filter_Biom")
kable(df)

	échantillon	mean	median	min	max
Biom	846	21339.08	5	0	2709474
Filter_Biom	846	21207.16	1092	181	2709474

#Info séquençage
#summary(Corres_ID$SeqDepth)
#median(Corres_ID$SeqDepth)

Le niveau d’abondance semble hétérogène. Il faudra donc prévoir un procédé de rarefaction

Graphique d’abondance des échantillons

• données complètes

plot_bar(Biom)+
  ggtitle("Abondance par échantillon sur données totales")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

• données filtrées

plot_bar(filter_biom)+
  ggtitle("Abondance par échantillon sur données filtrées")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

Dorénavant nous ne travaillerons que sur données filtrées Ne sert à rien d’aller au niveau Phylum ou infra car trop d’échantillons donc aucune couleur n’est visible Nous remarquons que le niveau d’abondance par échantillon est hétérogène de 25 000 en moyenne à 250 000.

Analyses effectuées mais retirées : J’avais fait des niveaux d’abondance par échantillon en fonction de différentes variables (sexe, satut tabagique…) mais ne sert à rien car ça équivaut à voir la profondeur et ce n’est pas biologiquement informatif donc enlevé.

Répartition des OTU par phylum

t <- table(tax_table(filter_biom)[,"Phylum"])
kable(t)

Var1	Freq
Actinobacteria	22
Bacteroidetes	223
Firmicutes	686
Fusobacteria	2
Proteobacteria	45
Synergistetes	9
Tenericutes	5
Verrucomicrobia	5

Visualisation de l’abondance en fonction des taxa

plot_bar(filter_biom, "Phylum", fill = "Phylum")

Les Firmicutes sont les bactéries majoritaires. Elles sont les plus abondantes suivies des Bacteroidetes

• Bar plot abondance

#Abondance par phylum
phy <- tax_glom(filter_biom, taxrank = "Phylum")
plotdata<-psmelt(phy)
ggplot(plotdata,aes(y=Abundance))+
  geom_boxplot(aes(group="Phylum", fill="Phylum"))+
  facet_wrap(~Phylum, scales="free_y",nrow=2)

• Abondance des firmicutes niveau famille au sein des échantillons

plot_composition(filter_biom, "Phylum", "Firmicutes", "Family", fill = "Family")+
  facet_grid(scales = "free_x") + 
  theme(axis.text.x = element_blank())

Pas de différence flagrante. Les échantillons semblent assez hommogène dans la composition par famille. Sur les données au niveau du genre rien ne se dégage ainsi qu’en fonction du sexe, age, tabac… Les échantillons semblent peu différents entre eux. Nous risquons donc d’avoir une faible voir une absence de diversité béta

Au vu des différents graphs, il y a une grande diversité de bactéries et nous ne voyons pas de genre majoritaire quelque soit la variable. Pour l’ordre Clostrdiales est majoritaires quelque soit la variable prise (exemple ci dessous)

Composition en Firmicutes au niveau de l’ordre en fonction de l’âge

plot_composition(filter_biom, "Phylum", "Firmicutes","Order",  fill = "Order") + 
  facet_grid(~AGE, scales = "free_x", space = "free_x")+
  theme(axis.text.x = element_blank())

Je pense que les graphs suivants ne sont d’aucune utilité . J’ai fait les mêmes sur les femmes mais trop d’échantillons

Alpha-diversité

Diversité intrinsèque à chaque échantillon

Etape de raréfaction

filter_rare <- rarefy_even_depth(filter_biom, rngseed = 20210416)
  
save(Biom, filter_biom, filter_rare, file = "Rarefied_Biom.Rdata")  

Graphique Covariable et test statistique

Exploration de l’impact de chaque covariable sur la diversité-alpha

Shannon, Cet indice représente à la fois le nombre d’espèces d’un milieu mais aussi la répartition des effectifs individuels au sein des espèces présentes. Simpson mesure de régularité cad mesure la probabilité que deux individus pris au hasard appartiennent à la même espèce Chao1 tient davantage compte des espèces peu abondantes nécéssite de conserver singleton donc pas pris

plot_richness(filter_rare, x = "AGE", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = AGE)) ## add one boxplot per AGE  values

# Couleur ne fonctionne pas scale_fill_manual(values=c("orange", "blue", "green3","red","yellow")) 

Age diffère significativement en terme d’observation d’OTU

div_data <- cbind(estimate_richness(filter_biom),  ## diversity indices
                  sample_data(filter_biom)         ## covariates
                  )

Mais également en terme de diversité en nombre d’espèce

plot_richness(filter_rare, x = "SEX", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = SEX)) ## add one boxplot per SEX  values

Pas de différence signicative pour le nbre d’OTU mais différence significative pour le nbre d’espèces

plot_richness(filter_rare, x = "tabac", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = tabac)) ## add one boxplot per AGE  values

Observation significative ainsi que Shannon. InvSimpson 0.08

Conclusion au vu des différents tests sur données filtrées non raréfiées:

Attention cependant dans notre cas d’alpha diversité sur la présence /abs (observed richness) il convient de travailler sur données raréfiées

• Pour le modèle Observé Tabac Age et BMI significatif (, ,)
• Pour l’indice de Shannon toutes les variables sont significatives voir hautement significatives (age SEX)
• Pour l’indice de InvSimpson les variables sauf tabac sont significatives voir hautement significatives (age SEX)
• Aucune interaction n’est significative

div_data <- cbind(estimate_richness(filter_rare),  ## diversity indices
                  sample_data(filter_rare))         ## covariates

model <- aov(Observed ~ tabac+age+SEX+BMI, data = div_data)
anova(model)

model <- aov(Shannon ~ tabac+age+SEX+BMI, data = div_data)
anova(model)

model <- aov(InvSimpson ~ tabac+age+SEX+BMI, data = div_data)
anova(model)

Conclusion au vu des différents tests sur données filtrées raréfiées:

• Pour le modèle Observé Tabac Age et BMI significatif (, ,) SEX ns
• Pour l’indice de Shannon toutes les variables sont significatives voir hautement significatives (age SEX)
• Pour l’indice de InvSimpson les variables sauf tabac sont significatives voir hautement significatives (age SEX)
• Aucune interaction n’est significative

Diversité Béta :

correspond à la différence de diversité des espèces entre plusieurs milieux Diversité entre échantillons : indice de dissimilarité (Bray et Curtis , de Jacard) Phylogénie (Unifrac) ### Nous n’avons pas d’arbre phylogénétique donc pas de distance unifrac et wunifrac

Ordination sans contrainte: capture de la diversité peu importe d’où elle provient.

Conclusion : Aucune tendance se dégage. Il n’y a pas de différence de diversité d’espèces entre les différent échantillons pour les catégories (age, sexe, tabac, bmi). Impossible de voir pour l’activité physique. La plage semble trop importante pour une coloration. Voir pour mettre des catégories

dist.jac <- distance(filter_rare, method = "jaccard")
dist.bc <- distance(filter_rare, method = "bray")
#dist.uf <- distance(xxx, method = "unifrac")
#dist.wuf <- distance(xxx, method = "wunifrac")

Représentation MDS (pour « Metric MultiDimensional Scaling » soit analyse multidimensionnelle métrique) La NMDS ne converge pas et NMDS déforme l’espace pour faire apparaitre des groupes éventuels. Du coup les distances apparentes ne sont pasfidèles aux distances réelles.

my_palette <- c('red','yellow','purple','green','orange')
p.jac <- plot_ordination(filter_rare, 
                         ordinate(filter_rare, method = "MDS", distance = dist.jac), 
                         color = "age") + 
  ggtitle("Jaccard")+ scale_color_manual(values = my_palette)



p.bc <- plot_ordination(filter_rare, 
                         ordinate(filter_rare, method = "MDS", distance = dist.bc), 
                         color = "age") + 
  ggtitle("Bray-Curtis")+ scale_color_manual(values = my_palette)

gridExtra::grid.arrange(p.jac, p.bc, ncol = 2)

Quelque soit la variable utilisée, aucun structuration n’est identifiée. Les 2 premiers axes ne capturent que que 22% de la diversité avec la distance de Bray-Curtis.

/ /

Analyses complémentaires non essentielles au vu des résultats précédents

Ordination avec contraintes: capture de la diversité issu de certaines covariables.

Risque de ne rien voir car déjà rien ne ressort sans contrainte Changement de l’argument method en CAP (Constrained Analysis of Proximities) et indiqué quelle covariable peut expliquer cette diversité

ord <- ordinate(filter_rare, method = "CAP", distance = dist.bc, formula = ~ age)
p12 <- plot_ordination(filter_rare, ord, color = "age", axes = c(1:2))
p23 <- plot_ordination(filter_rare, ord, color = "age", axes = c(2:3))
gridExtra::grid.arrange(p12, p23, nrow = 1)

/

Clustering Hierarchique

Aucune clusterisation des échantillons n’est visible Si présence d’arbre phylo possibilité de faire avec en utilisant dist=“unifrac”

par(mar = c(1, 0, 2, 0))
plot_clust(filter_rare, dist = "bray", method = "ward.D2", color = "SEX", 
           title = "Clustering of samples (Bray-Curtis + Ward)\nsamples colored by sexe")

/ /

PERMANOVA

metadata <- sample_data(filter_rare) %>% as("data.frame")
model <- vegan::adonis(dist.bc ~ age + SEX + tabac + BMI + APhysGlobHParSem , data = metadata, permutations = 999)
model

Call:
vegan::adonis(formula = dist.bc ~ age + SEX + tabac + BMI + APhysGlobHParSem,      data = metadata, permutations = 999) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                  Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
age                4     1.732 0.43289  2.0403 0.00954  0.001 ***
SEX                1     1.273 1.27315  6.0007 0.00701  0.001 ***
tabac              2     0.672 0.33604  1.5839 0.00370  0.013 *  
BMI                1     0.318 0.31767  1.4973 0.00175  0.065 .  
APhysGlobHParSem   1     0.198 0.19769  0.9318 0.00109  0.552    
Residuals        836   177.370 0.21216         0.97691           
Total            845   181.562                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Toutes les variables, à l’exception de l’activité physique, sont significatives. Cependant elles n’expliquent quasiment pas de variabilités < 1%

/ / ## Heatmap

Ne voyant pas de cluster, il n’est pas nécéssaire de faire une heatmap d’une variable en considérant les distances de Bray-Curtis.

plot_heatmap(filter_rare, low = "yellow", high = "red", na.value = "white") + 
  facet_grid(~age, scales = "free_x", space = "free_x")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

Quelque soit la tranche d’âge, il y a beaucoup de diversité. Ce sont les mêmes OTU qui sont présents dans tous les échantillons quelque soit la catégorie d’âge.

 # Possibilité de zoomer sur une zone les 50 taxa les plus élevés

top_50_taxa <- taxa_sums(filter_rare) %>% sort(decreasing = TRUE) %>% names() %>% head(n = 50)
p <- plot_heatmap(prune_taxa(top_50_taxa, filter_rare), 
             low = "yellow", high = "red", na.value = "white") + 
  facet_grid(~age, scales = "free_x", space = "free_x")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
  
plot(p)

Anlayses différentielles

Objectif : Voir si l’abondance de certaines espèces diffèrent entre les groupes

knitr::opts_chunk$set(fig.height = 5)

cds <- phyloseq_to_deseq2(filter_biom, ~SEX)
dds <- DESeq2::DESeq(cds, sfType = "poscounts")   
results <- DESeq2::results(dds) %>% as.data.frame()

##### On ne peut prendre que les OTU dont le fold change >0.05% on peut être plus drastique
da_otus <- results %>% as_tibble(rownames = "OTU") %>% 
  filter(padj < 0.05) %>% 
  arrange(log2FoldChange) %>% 
  pull(OTU)
length(da_otus)

[1] 60

##### Récupères ces OTU significatifs et on les ordonne dans une heatmap
p <- plot_heatmap(prune_taxa(da_otus, filter_biom), ## keep only da otus...
                  taxa.order = da_otus       ## ordered according to fold-change
                  ) + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ 
  facet_grid(~SEX, scales = "free_x", space = "free_x") + 
  scale_fill_gradient2(low = "#1a9850", mid = "#ffffbf", high = "#d73027",
                       na.value = "white", trans = log_trans(4),
                       midpoint = log(100, base = 4))

plot(p)

cds <- phyloseq_to_deseq2(filter_biom, ~age)
dds <- DESeq2::DESeq(cds, sfType = "poscounts")   
results <- DESeq2::results(dds) %>% as.data.frame()

##### On ne peut prendre que les OTU dont le fold change >0.05% on peut être plus drastique
da_otus <- results %>% as_tibble(rownames = "OTU") %>% 
  filter(padj < 0.05) %>% 
  arrange(log2FoldChange) %>% 
  pull(OTU)
length(da_otus)

[1] 57

##### Récupères ces OTU significatifs et on les ordonne dans une heatmap
p <- plot_heatmap(prune_taxa(da_otus, filter_biom), ## keep only da otus...
                  taxa.order = da_otus       ## ordered according to fold-change
                  ) + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  facet_grid(~age, scales = "free_x", space = "free_x") + 
  scale_fill_gradient2(low = "#1a9850", mid = "#ffffbf", high = "#d73027",
                       na.value = "white", trans = log_trans(4),
                       midpoint = log(100, base = 4))



plotly::ggplotly(p)

Conclusion

• Beaucoup de bactéries dominée principalement poar Firmicutes et Bacteroidetes.
  • Pas de genre spécifique à une catégorie
  • Grande diversité de bactéries au sein de chaque échantillon. Différences significatives pour le nombre d’espèces de bactéries entre les personnes de différentes catégories Les différentes variables sont significativement différentes
  • Homogéniété en terme de représentativité des bactéries au sein des échantillons